%  This file is part of Erlywog.

%  Erlywog is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.

%  Erlywog is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.

%  You should have received a copy of the GNU General Public License
%  along with Erlywog.  If not, see <http://www.gnu.org/licenses/>.

%  Copyright Jasm Sison 2009

-module(mersenne).

-export([new/1, generate/1, unbox/1, rand/2]).

%% The coefficients for MT19937 are:
-define(W, 32). %% word size (in bits)
-define(N, 624). %% degree of reccurence
-define(M, 397). %% middle word, 1 <= M <= N
-define(R, 31). %% seperation of word, 0 <= R <= W-1
-define(A, 16#9908B0DF). %% coefficients of the rational normal form twist matrix
-define(B, 16#9D2C5680). %% TGFSR(R) tempering bitmasks
-define(C, 16#EFC60000). %% TGFSR(R) tempering bitmasks
-define(S, 7).  %% TGFSR(R) tempering bit shifts
-define(T, 15). %% TGFSR(R) tempering bit shifts
-define(U, 11). %% additional Mersenne Twister tempering bit shifts
-define(L, 18). %% additional Mersenne Twister tempering bit shifts
-define(INIT_NUMBER_MAGIC, 1812433253). %% determined by the algorithm inventor, do not change!
-define(INIT_RIGHT_SHIFT, 30). %% another magic number

%% accepts a random seed of any value > 0, recommands a prime number
%% returns MT = {index = 0, mersenne array}
new(Seed) ->
   A = lists:foldl(fun (0, InitVal) ->
                          [<<InitVal:?W>>];
                       (I, Acc = [<<O:?W>>|_T]) ->
                          [<<(?INIT_NUMBER_MAGIC * (O bxor (O bsr ?INIT_RIGHT_SHIFT)) + I):?W>>|Acc]
               end, Seed, lists:seq(0,?N-1)),
   {0, lists:reverse(A)}.

%% generate N random integers
rand(N, MT) ->
   [H|NN] = lists:foldl(fun (_I, [MT1|Ns]) ->
                  {R, MT2} = generate(MT1),
                  L = [unbox(R)|Ns],
                  [MT2|L]
               end, [MT], lists:seq(1, N)),
   {lists:reverse(NN), H}.


twist(<<Y0:?W>>) ->
   Y1 = Y0 bxor (Y0 bsr ?U), %%shift 11
   Y2 = Y1 bxor ((Y1 bsl ?S) band ?B), %% shift by 7 and 0x9d2c5680
   Y3 = Y2 bxor ((Y2 bsl ?T) band ?C), %% shift by 15 and 0xefc60000
   <<(Y3 bxor (Y3 bsr ?L)):?W>>. %% shift by 18
     

%% accepts an MT generate by mersenne:new/1
%% returns MT = {index = 0, mersenne array}
generate({0, A}) ->
   MT = step(A),
   {twist(hd(MT)), {1 rem ?N, MT}};
generate({I, A}) ->
   {twist(lists:nth(I+1, A)), {(I + 1) rem ?N, A}}.

lower_bit(<<_:1, N:(?W - 1)>>) -> N.
   
upper_bit(<<N:1, _:(?W - 1)>>) -> N.

unbox(<<N:?W>>) -> N.
               
step(MTL) ->
   MT1 = array:from_list(MTL),
   MT2 = array:foldl(fun (I, _Val, MT) ->
                           Y = upper_bit(array:get(I, MT)) + lower_bit(array:get((I+1) rem ?N, MT)),
                           Y1 = unbox(array:get((I+?M) rem ?N, MT)) bxor (Y bsr 1),
                           Y2 = if (Y rem 2) == 1 -> %% is y odd?
                                      Y1 bxor ?A;
                                   true ->
                                      Y1
                                end,
                           array:set(I, <<Y2:?W>>, MT)
                     end, MT1, MT1),
   array:to_list(MT2).

